Thermodynamic properties of the two-dimensional 5 = 1/2 Heisenberg antiferromagnet 

coupled to bond phonons 
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By applying a quantum Monte Carlo procedure based on the loop algorithm we investigate ther- 
modynamic properties of the two-dimensional antiferromagnetic 5=1/2 Heisenberg model coupled 
to Einstein phonons on the bonds. The temperature dependence of the magnetic susceptibility, 
mean phonon occupation numbers and the specific heat are discussed in detail. We study the spin 
correlation function both in the regime of weak and strong spin phonon coupling (coupling constants 
g = 0.1, U) = 8 J and g — 2, u) = 2 J, respectively). A finite size scaling analysis of the correlation 
length indicates that in both cases long range Neel order is established in the ground state. 
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I. INTRODUCTION 

During the last years, there has been considerable in- 
terest in low-dimensional spin systems with spin phonon 
coupling. While one-dimensional (ID) models including 
this mechanism have been studied extensively, little is 
known about two-dimensional (2D) systems with spin 
phonon coupling. 

In ID systems, the mechanism of spin phonon cou- 
pling is closely connected with the phenomenon of the 
spin Peierls transition. Theoretical understanding of this 
phase transition goes back to the work of Pytte^i who 
showed that a three-dimensional system consisting of 
uniform antiferromagnetic S =1/2 Heisenberg chains is 
unstable towards dimerization of the chains if coupled 
to three-dimensional lattice vibrations. Additionally, he 
showed that in the adiabatic limit of small phonon fre- 
quencies the dimerized phase can be described by a stat- 
ically dimerized spin model with temperature-dependent 
dimerization. Cross and Fished improved on the calcu- 
lation of Ref. Q] by treating the spin part of the Hamilto- 
nian in continuum field theory* Their calculation yielded 
a convincing description of the phonon softening in the 
limit of small spin phonon couplings. In the case of 
CuGeOa, however, which was the first inorganic spin 
Peierls compound discovered^ this treatment is not suf- 
ficient. For this reason various ID dynamical models 
have been investigated, taking into account the coupling 
of magnetic and phononic degrees of freedom explicitly. 
Most publications on this issue refer to two different 
models, both showing a quantum phase transition be- 
tween a dimerized and a Neel ordered phase. While in 
the so-called difference coupling mode&Si2iS the mag- 
netic interaction between nearest neighbours depends on 
the distance between neighboured sites, the bond cou- 
pling model£ii2*ii*i2iiii4 is considered to be more realis- 
tic for describing the spin phonon coupling mechanism in 
CuGe0 3 »i^ 

In the case of the 2D 5 = 1/2 Heisenberg antiferromag- 
net with spin phonon coupling discussions of statically 
dimerized models are found in the literaturepiiiSiiLSLSi 
As in one dimension, these models are thought to de- 



scribe the dimerized phase of dynamical models in the 
adiabatic limit. In contrast to the ID case though it is 
not clear how to place alternating magnetic couplings in 
both spatial directions on a square lattice. By comparing 
ground state energies of three statically dimerized models 
Sirker et. alm^ conclude that a stair-like distortion of the 
lattice is the energetically favoured dimerization pattern, 
contradicting an older result by Tang and Hirschii who 
find a plaquette-like distortion. 

The aim of this paper is to investigate thermodynamic 
properties of the 2D Heisenberg model coupled to bond 
phonons. Such a model takes into account the elastic 
energy due to lattice distortions which is not included 
in statically dimerized models. Furthermore, the whole 
range of phonon frequencies is accessible in our treat- 
ment. 

The paper is organized as follows. In Sect. [H]we in- 
troduce the model Hamiltonian and describe our quan- 
tum Monte Carlo algorithm. In Sect. IIIII we discuss the 
temperature dependence of the magnetic susceptibility, 
mean phonon occupation numbers and the specific heat. 
An analysis of the spin correlation function is found in 
Sect. IIVI Section fVl concludes with a summary. 



II. MODEL HAMILTONIAN AND QUANTUM 
MONTE CARLO METHOD 



The model we consider is a genera lization of the ID 
bond coupling model from R.efs. Illlldl The Hamiltonian 
readsS^ 
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Here <xy denote Pauli spin operators at lattice site 
on a square lattice, while ay and ay- (by and by) are 
phonon annihilation and creation operators on the bond 
between site and site (i + (between sites 
and (i,j + 1)). Note that we assume periodic boundary 
conditions. 

By shifting the phonon operators according to 
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and neglecting a constant energy contribution model JJJ 
can be mapped onto the phenomenologically more real- 
istic Hamiltonian 

H = - ^2(J' + g'[aij + al^SijBi+ij 

ij 

+ I Y,( J> + 3%j + ( 3 ) 

ij 

+ ^^2( a ij a ij + b ij°ij)> 
ij 

with rescaled coupling constants J' = J(l + g 2 J/oj) 
and g' = gj. Note that Eq. J3J differs from 
JTJ by the absence of the (unphysical) static terms 
-Jg/2(alj + a lj )-Jg/2(b\ J + by). Model ©, however, 
cannot be analyzed directly in a quantum Monte Carlo 
study because a sign problem occurs. For this reason we 
choose Hamiltonian P as a starting point for our anal- 
ysis and shift the numerical results if needed. 

To study the properties of model (J), we developed a 
quantum Monte CarloS^ algorithm similar to the algo- 
rithm described in Ref. 0. First, the partition function 
of the 2D quantum system is mapped onto a three- 
dimensional classical system. Technically, this is done by 
means of a Trotter Suzuki decomposition^! We then ap- 
ply an update procedure which treats spin and phononic 
degrees of freedom separately. For the spin updates, we 
make use of a modified loop algorithm2a2^ for quantum 
spin systems. The main advantage of the loop algorithm 
is that it allows global spin updates, substantially re- 
ducing autocorrelation times. Furthermore, so-called im- 
proved estimators can be used in the evaluation of the 
magnetic susceptibility and spin correlations. To mod- 
ify the phonon occupation numbers we apply local heat 
bath updates. By building clusters of phonons in imagi- 
nary time direction we extended the algorithm from Ref. 
ITU diminishing autocorrelation effects even more. Obvi- 
ously the detailed balance condition is fulfilled for both 
steps separately and thus for the whole procedure. 

For the phonon updates we had to introduce a cut- 
off, allowing occupation numbers up to 40 phonons per 
bond. The effect of such a truncation of the Hilbert space 
is negligible if the measured mean phonon occupation 
numbers are more than an order of magnitude smaller 
than the cutoff. In Sect. IIII Bl we show explicitly that 
this condition is fulfilled. In order to take into account 
the high dimension of the phonon subspace of the Hilbert 



space, we employed the importance sampling technique 
and made 30 phonon updates per spin update, using only 
the last configuration for the evaluation of expectation 
values. Additionally, for each temperature the first 25% 
of the sweeps were skipped for thcrmalization. 

For Monte Carlo simulations based on a Trotter Suzuki 
decomposition the estimates of thermodynamic quanti- 
ties depend on the inverse Trotter number squared^ In 
the following sections we give the explicit value for the 
Trotter number M. With the values for M chosen, we 
find the statistical fluctuations of our results larger than 
the effect of the finite Trotter number. 

Before discussing the results in detail we add one fur- 
ther remark concerning statistical errors. In our cal- 
culations we neglected autocorrelation effects and - as 
a rough estimate - calculated root-mean-squared errors 
only. If no error bars in the plots of this paper are shown, 
the errors are smaller than the symbol size used. 



III. THERMODYNAMIC PROPERTIES 

In this section we discuss the finite temperature prop- 
erties of model JQ|. We expect that the knowledge of 
how a non- vanishing spin phonon coupling influences the 
thermodynamic properties will be of importance for the 
interpretation of experiments. This might be of partic- 
ular interest for substances which display e. g. acoustic 
anomalies or for which it is known that the exchange 
integral depends sensitively on the positions of the ions. 

Here, we confine ourselves to temperatures 
0.5 J < T < AJ . In this temperature range, we find 
that the dependence of measured quantities on the 
system size is negligible if we consider linear system 
sizes N > 12. All results presented in this section were 
calculated on a lattice with 12x12 sites, providing state- 
ments about system properties in the thermodynamic 
limit. If not stated differently, at each temperature 10 5 
spin updates were executed. For the Trotter number a 
value of M = 80 was chosen. 



A. Magnetic susceptibility 

We start with a discussion of the magnetic suscepti- 
bility per site x f° r vanishing magnetic fields. Figure ^ 
shows the dependence of the susceptibility on the spin 
phonon coupling g for fixed phonon frequency ui, Fig. [21 
the dependence on uj for fixed g. In both figures Monte 
Carlo results for the susceptibility of the 2D Heisenberg 
model are included. 

The results can be summarized as follows. For fixed ui, 
the overall height of the susceptibility is diminished with 
increasing spin phonon coupling. As in the ID case^i a 
large spin phonon coupling tends to reduce the magnetic 
response of the system. On the other hand, for fixed 
g the susceptibility is growing with increasing phonon 
frequency. In the antiadiabatic limit phononic degrees 
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FIG. 1: Magnetic susceptibility vs. temperature for fixed 
ui — 2J and spin phonon coupling g between 0.5 and 2.0. For 
comparison Monte Carlo results for the Heisenberg model are 
plotted. 
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FIG. 2: Magnetic susceptibility vs. temperature for fixed 
g = 1 and phonon frequencies between 2J and 8J. Again 
for comparison Monte Carlo results for the Heisenberg model 
are shown. 



FIG. 3: The data from Fig.Qin terms of the rescaled magnetic 
coupling J'. 



0.08 



0.06 



0.04 




0.02 /»» 



— ID Heisenberg model, exact solution _ 
o-o ID, g=0.5, C0=2J 
b--h 2D Heisenberg model 
o-o 2D, g=0.5, m=2J 



2 

T/J 



FIG. 4: Temperature dependence of the magnetic susceptibil- 
ity in d = 1 and d = 2 for gr = 0.5 an d a) = 2 J. For comparison 
in d = 1 the exact result from Refs. l27l2Sl and in d = 2 Monte 
Carlo data for the Heisenberg model are shown. 



of freedom are suppressed, yielding the results of the 
Heisenberg model. The curves show a broad maximum 
which is typical for antiferromagnetic spin models. This 
maximum is shifted to higher temperatures with increas- 
ing gj/u. 

We find that both the shift of the position of the max- 
imum and the reduction of magnetic response with in- 
creasing is mainly due to the static terms in (JIJ. In 
units of the rescaled magnetic coupling J' from the trans- 
formed Hamiltonian @ the shift of the maximum is re- 
duced significantly (see Fig.[3J), and the reduction of the 
magnetic response due to the spin phonon coupling is 
not very strong. The same behaviour has already been 
reported for the ID bond coupling model in Ref. Il4l 



For comparison between d = 1 and d — 2 we return 
to model without phonon shift. Figure 01 shows the 
susceptibilities for g = 0.5, lo = 2J. Compared to the 
ID case, the overall height of the susceptibility of the 
2D model is diminished. This effect is explained by the 
larger coordination number in the 2D case, reducing the 
response of the system to an external magnetic field. 

Qualitatively the influence of the spin phonon coupling 
is similar in d = 1 and d = 2. In Fig. 0] the exact result 
from Refs. for the ID and the Monte Carlo results 

for the 2D Heisenberg model are shown. Both in ID 
and 2D we find a significant shift of the maximum and a 
strong reduction of the maximum height as compared to 
the Heisenberg model. 
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B. Mean phonon occupation numbers 

Further insight into the model can be gained by study- 
ing the influence of the spin phonon coupling on the mean 
phonon occupation numbers 



2.5 



(n) 
(m) 



ij 



(4) 
(5) 



as compared to the free phonon case. These numbers 
can be viewed as a measure of the strength of lattice 
vibrations and therefore allow to analyze how the lattice 
is influenced by the spin degrees of freedom. 

As expected, we find no difference in the mean occu- 
pation numbers (n) and (m), and therefore restrict the 
following discussion to the behaviour of (n) . Figure |S] 
shows the Monte Carlo results for (n) for different values 
of g and w in a plot vs. T/ui. The data are compared to 
the Bose distribution for free Einstein phonons 
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which is also shown in Fig. Again we find a striking 
similarity to the ID bond coupling model. In d = 1 it has 
been found that the mean phonon occupation numbers 
obey the relation™. 



(n)(T) =n + nfree(T) 



(7) 



with a temperature-independent constant uq. As a good 
approximation, this relation is valid in a temperature 
range 0.5 J < T < 3 J in the 2D case as well. This can be 
seen most clearly in Fig. [SI which shows the same results 
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FIG. 6: Differences (n) — nf rcc vs. T/u 



as Fig. [3] with relation @ subtracted from the Monte 
Carlo data. 

In order to derive an expression for n , we averaged 
the differences (n) — n{ TCC over the temperature range 
0.5J < T < 3 J and plotted these values vs. g 2 J 2 /uj 2 (see 
Fig. [7J|. By applying linear regression we find that in 
d = 2 the shift obeys the relation 



n w (1.375 ±0.003) ( — 



(8) 



Thus the only difference between the relations for the 
mean phonon occupation numbers in the ID and 2D case 
is given by the numerical prefactor in (jSJ) , the value being 
1.375 in d = 2 and 2 in d = l»±i 

We close this section with a technical remark concern- 
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FIG. 5: Mean phonon occupation numbers (n) as a function 
of T/uj. The solid line shows the Bose distribution rif ree for 
free Einstein phonons. 



FIG. 7: Over the temperature range 0.5 J < T < 3J averaged 
values no vs. g 2 J 2 /uj 2 . The solid line shows the result from 
linear regression. 
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ing our choice for the cutoff for the phonon occupation 
numbers. As can be seen in Fig. [21 in the whole temper- 
ature range the measured mean occupation numbers are 
more than an order of magnitude smaller than the cutoff 
40. In retrospect our choice therefore is justified. With 
(JZJ) and (JHJ) we have also found an expression that might 
be important for other numerical methods (e. g. exact di- 
agonalization) which depend crucially on a (low) cutoff 
in the phonon numbers. 



C. Specific heat 

Another important thermodynamic quantity is the 
specific heat per site C. Although in principle this ob- 
servable is directly accessible in the experiment, it is dom- 
inated by lattice vibrations, making it difficult to extract 
its magnetic part. Even for a simple model as given by 
Hamiltonian Q, we find this behaviour confirmed. Fig- 
ure [SI shows Monte Carlo data for the specific heat in a 
system with g — 2 and uj — 2 J and the exact result 

C frcc (T) = 2 (-) (e , /T _ 1)2 (9) 

for free phonons of the same frequency (the factor two 
accounts for two phonons per lattice site). There is only 
a small difference in the overall height of C and Cf ree . 
At high temperatures, both curves approach the same 
constant value, yielding the Dulong-Petit rule. 

The spin phonon coupling influences the specific heat 
significantly though. As can be seen in Fig. [5] as well, the 
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FIG. 8: Specific heat plotted vs. temperature for g = 2 and 
uj = 2J. The solid line shows result @ for free phonons of 
the same frequency. Plotted are also Monte Carlo results for 
the 2D Heisenberg model and the sum of the free phonon and 
Heisenberg results. The number of spin updates is 5 x 10 5 for 
the system with spin phonon coupling and 10 5 in the Heisen- 
berg case. 



curve for the system with spin phonon coupling differs 
significantly from the sum of Monte Carlo results for the 
2D Heisenberg model and the contribution © for free 
Einstein phonons. Note that both the strong fluctuations 
and the divergency of the data for T — > are due to 
difficulties in evaluating the specific heat within Monte 
Carlo procedures as discussed in Ref. |29j 

IV. SPIN CORRELATION FUNCTION AND 
GROUND STATE PROPERTIES 

We now turn our attention towards ground state prop- 
erties of model ||TJ . In principle the Monte Carlo method 
is only applicable at finite temperatures. By analyz- 
ing the behaviour of the spin correlation function at low 
temperatures, however, it is possible to make statements 
about system properties at T = 0. 

The argument is as follows. Suggest that we choose 
the coupling constants in (|TJ such that the system is 
Heisenberg-like, showing long range Neel order in the 
ground state. Then for reasons of universality we expect 
that the spin correlation function 

r 

obeys the result known for the 2D Heisenberg 

G(d) ~ (-l) d i+^ \d\- x e-^/^), (11) 

with the algebraic exponent A closed! to the classical 
Ornstein-Zernike value of \. Here £(T) is the spin corre- 
lation length which can be interpreted as the mean size 
of domains with antiferromagnetic order. At T = these 
domains get macroscopic, because for T — » the corre- 
lation length diverges exponentially. 2£L2ii22i22i Assuming 
A = \ in 111(1 , this means that the static structure factor 

S{q) = Y f ^G{d) (12) 

d 

diverges for momentum 45 q = (tt, it) with the linear sys- 
tem size like A 2 . As long as the correlation length in the 
infinite system stays significantly larger than the system 
sizes considered, this behaviour should be visible at low 
temperatures. 

We first illustrate this in case of the 2D Heisenberg 
model (g = ui = in JTJ) at T = 0.1 J, taking 10 5 spin 
updates and choosing a Trotter number of M = 160 for 
the evaluation of spin correlations. At this temperature, 
the correlation length in the infinite system is of the order 
of 10 9 lattice spacings^S^iiSS^ As can be seen in the 
upper plot of Fig. the static structure factor shows a 
pronounced maximum at q = (w, w), and the peak height 
roughly scales with the system size like N*. Note that 
deviations from the A^-dependence might indicate that 
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FIG. 9: Static structure factor S(q, q) of spin correlations as 
a function of momentum q for a diagonal cut through the 
first Brillouin zone at T = 0.1J' for different system sizes. 
Shown are results for the 2D Heisenberg model (top) and for 
model Q with g = 0.1, lu = 8J (middle) and g = 2, lu — 2J 
(bottom), respectively. In all plots the inset shows the height 
of the maximum at (tt, tt) as a function of JV 5 . 



in the quantum system A differs slightly from the value 
of \ chosen above. 

We now return to model with spin phonon cou- 
pling and discuss our results for two different choices 
of coupling constants. To compare our results for dif- 
ferent values of g and u) we drop the unphysical terms 
-Jg/2(a\ j + dij), -Jg/2(b\ J + 6 y ) in Q. We therefore 
measure the temperature in units of the rescaled mag- 
netic coupling J' of the effective Hamiltonian J2J. We 
calculated spin correlations at T = 0.1J', taking 5 x 10 s 
spin updates and M — 160 in our calculations. First, we 
consider a system with a small value for ^- (g = 0.1, 
lu = 8J) where no dimerization is expected. Again 
the static structure factor shows a pronounced peak for 
q = (tt, 7r), and the peak height scales with the system 
size like TV 2 (middle of Fig. indicating Heisenberg-like 
behaviour as anticipated. 

For the second system we choose g = 2 and lu — 2J. 
In the ID case this choice corresponds to a system which 
strongly dimerizes in the ground stated Even here we 
find a pronounced peak of S for (tt, tt) (bottom of Fig.[5J), 
the maximum height scaling like The interpretation 
is that even in the regime of large values for the system 
shows antiferromagnetic order in the ground state. This 
is a striking difference to the ID bond coupling model. 
Compared to the case with g = 0.1 and lu — 8 J, however, 
we find that the spin phonon coupling counteracts the 
tendency of the system to order antiferromagnetically in 
the ground state. As can be seen in Fig. for fixed 
system size the peak heights S(tt, tt) in the case of strong 
spin phonon coupling arc slightly diminished as compared 
to the weak coupling regime. 

Our results can be confirmed by a direct analysis of the 
temperature dependence of the spin correlation length. 
For both systems and at various temperatures we ex- 
tracted finite system correlation lengths £n by fitting the 
function 



/(d) = a(-l) 



di-\-d.2 



-(AH <*!)/« 




(13) 



with two free parameters a, £ to our data for system sizes 
N = 10,12,14,20,24. For g = 0.1 and w = 8 J, we se- 
lected values M = 120 for temperatures 0.5 J < T < 0.9 J 
and M = 80 for T > J, taking 10 5 spin updates (for 
N = 24 we chose M = 120 for all temperatures). For 
9 = 2 and ui — 2J the choice was M — 160 for 
1.5J < T < 1.9J, M = 120 for 2J < T < 4J and 
M — 80 for T > 5J, again averaged over 10 5 Monte 
Carlo sweeps (for N = 20 we took 1.5 x 10 5 , for N = 24 
and 1.5J < T < 1.9J we took 2 x 10 5 spin updates). 
The values for the Trotter number are sufficiently large 
to avoid finite size effects in Trotter direction. 

As has been said above, in case of the Heisen- 
berg model in leading order £ behaves like eT at low 
temperaturesi^2i2i*SiS In the upper panel of Fig. 1101 the 
natural logarithm of £at is plotted vs. y for g = 0.1 and 
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FIG. 10: Natural logarithm of finite system correlation 
lengths £jv vs. inverse temperature for five different system 
sizes for g = 0.1, uj = 8J (top) and g — 2, uj — 2J (bottom), 
respectively. Note that all temperatures are given in units of 
J'. 



lu = 8J. At high temperatures, no dependence on the 
system size is visible, and as expected we find the same 
behaviour as in the Heisenberg-model. At low tempera- 
tures finite size effects become important and the curves 
branch off from the asymptotic linear behaviour. In the 
case of strong spin phonon coupling the graph shows very 
similar features (bottom of Fig. HOfl . We therefore find 
Heisenberg-like behaviour of £ at finite temperatures also 
in the regime of strong spin phonon coupling. The main 
difference between the two cases is that at the same effec- 
tive temperature the correlation length for g = 2, uj = 2 J 
is significantly smaller than for g — 0.1, uj — 8 J. The 
analysis of £ therefore also implies that a strong spin 
phonon coupling weakens antiferromagnetic order. Both 
observations are consistent with the conclusions drawn 
from the analysis of S(q). Note that in both plots of 
Fig- El the temperatures are given in units of J'. 

By means of scaling arguments our analysis can be 
extended to make direct statements about ground state 



FIG. 11: Test of the scaling prediction l|14|l for g — 0.1, 
uj — 8 J (top) and g — 2, uj — 2 J (bottom). 



properties. Suppose the system shows long range Neel or- 
der in the ground state. In terms of the renormalization 
group this means that there is a critical fixed point at 
T = which controls the system properties at low tem- 
peratures. In this case a finite size scaling ansat32&2i2£ 



= F 



N 



(14) 



holds, where J 7, is a universal scaling function. With the 
data from Fig. ^| it is possible to test the scaling pre- 
diction JT3J. Plotting ^ vs. with N = 10,12 we 
find that both for g = 0.1, u> = 8 J (top of Fig. [nj and 
g = 2, uj = 2 J (bottom of Fig. 111(1 the data lie on one 
curve. The shape of the scaling function F in l|14[) . how- 
ever, depends on the choice of the coupling constants. 
The interpretation is that in the weak and in the strong 
coupling regime model (JIJ shows long range Neel order 
in the ground state, strongly confirming the conclusions 
drawn from our analysis of the static structure factor at 
low temperatures. 
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FIG. 12: Static structure factors S 1Y> (q) of spin correlations 
vs. momentum q in the ID case for different system sizes and 
temperature T = 0.1J'. In the upper graph the coupling 
constants are g = 0.1 and u) = 8 J, in the lower graph g — 2 
and uj = 2J. In both plots the inset shows the height of the 
maximum S 1D (ir) as a function of TV. 

It is instructive to compare these results to the low 
temperature behaviour of spin correlations in d = 1. Our 
argumentation is completely analogous to the 2D case. 
As has been discussed in Sect. the ID model shows 
a quantum phase transition between a Neel ordered and 
a dimerized phaseAii*^ Though only the approximate 
shape of the phase separation line in coupling constant 
space has been determined, it is known that for small 
values — the system is Heisenberg-like, showing quasi 
long range Neel order in the ground state. This means 
that the ID correlation function 

1 N 

G m {d) = -Y,(Wi+d) (15) 

i=i 

decays exponentially at finite temperatures, with a rate 
given by the correlation length £ 1D oc i. ,36,37,38,39,40 At 
T = 0, there is a crossover to an algebraic decay 3 ^*^^ 



G (d) ~ (16) 
and at T — the static structure factor 

N 

S 1B (q) =Y,e iqd G m (d) (17) 

diverges for momentum q = tt like the 7Vth partial sum 
of the harmonic series with the system size. As in d = 2, 
signs of this divergence should be visible at low temper- 
atures. 

For large values — , on the other hand, the chain 
dimerizes. This means that long range dimer order is 
established in the ground state. Therefore the spin corre- 
lation length £ 1D stays finite at T — 0, and we expect no 
dependence of 5* 1D on N in the low temperature regime. 

In Fig.^Jthe static structure factors S 1D for two sys- 
tems with the same choice of coupling constants as in 
d = 2 are plotted. As in d = 2, we measure the tem- 
perature in units of the rescaled magnetic coupling J' of 
the ID counterpart of J3J. We selected T = 0.1 J' and 
M = 160, executing 5 x 10 5 spin updates. In both sys- 
tems S 1T> shows a maximum at q = tt. In the Heisenberg- 
like system the maximum is more pronounced though, 
and for small system sizes the peak height depends on 
N. For larger system sizes, however, such a behaviour is 
not visible. This leads to the conclusion that in contrast 
to the 2D case the correlation length in the infinite sys- 
tem is not larger than the system sizes in consideration. 
In the system with dimerization in the ground state we 
find the expected behaviour: The values S 1d (tt) do not 
depend on the system size, indicating that £ 1D is very 
small. 

By analyzing the temperature dependence of £ 1D it 
is possible to distinguish more clearly between the two 
regimes. Similar to d = 2 we extracted correlation 
lengths £ 1D by fitting an exponential decay with two free 
parameters to our data. For g = 0.1 and u — 8 J, we 
executed 5 x 10 5 spin updates and selected M = 160, 
which is large enough to avoid effects by the finite Trot- 
ter number. Furthermore, the system sizes were chosen 
that large that finite size effects are negligible (N = 500 
for T = 0.025 J, 0.05 J, N = 400 for T = 0.075 J, 0.125 J 
and N = 300 for T = 0.1J, 0.15J). This can be seen 
in the inset of Fig. ^3 where is plotted vs. N at 
T = 0.05J w 0.05J'. For g = 2 and w = 2J also 5 x 10 5 
spin updates were made. The correlation lengths are that 
small that a chain length of N = 200 is sufficient to make 
statements about the thermodynamic limit. The effect of 
the finite Trotter number is more important in this sys- 
tem. For M — 400 at the lowest temperatures, however, 
the effect is smaller than the error which enters our anal- 
ysis during the fitting procedure (we chose M — 400 for 
T = 0.05J,0.075J,0.1J, M = 360 for T = 0.15 J and 
M = 160 for T > 0.2 J). 
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FIG. 13: Correlation lengths £ 1D vs. inverse temperature in 
the ID case. Note that the temperatures are given in units 
of J'. The inset shows the finite size behaviour of £jp for 
g = 0.1, lj = 8J at T = 0.05J'. 



Figurc^Dshows the correlation lengths in a plot vs. the 
inverse temperature in units of the rescaled effective cou- 
pling J'. In both systems we find the expected be- 
haviour: In the Heisenberg-like system the correlation 
length grows linearly with the inverse temperature, while 
in the dimerized system £ 1D takes small values and shows 
no such dependence in the temperature range shown. 

We close this section with a final remark concerning 
the 2D model JTJ. The results from this section need 
not mean that the model shows no lattice distortion. In 
both the statically dimerized stair and plaquette models 
e. g. a phase with coexisting dimerization and long range 
antiferromagnetic order is known^i For model Q it is 
therefore conceivable that small lattice distortions appear 
which cannot be detected by analyzing spin correlations 
at low temperatures. Even a finite temperature phase 
transition seems possible, because due to the Mermin- 



Wagner theorem^ a breaking of the discrete lattice sym- 
metry at finite temperatures cannot be excluded in a 2D 
system. However, we find no hints on a finite tempera- 
ture phase transition in the behaviour of thermodynamic 
properties in the temperature range discussed in Sect. 
11111 Further investigations of the model seem appropri- 
ate to clarify whether dimerization appears and which 
dimerization pattern is realized. 



V. SUMMARY 

By combining loop updates for spin and cluster up- 
dates for phononic degrees of freedom we have developed 
a quantum Monte Carlo algorithm to study the proper- 
ties of the 2D antiferromagnetic Heisenberg model cou- 
pled to bond phonons. 

As thermodynamic quantities are concerned, we 
studied the susceptibility, mean phonon occupation 
numbers and specific heat in the temperature range 
0.5 J < T < 4J. The properties of the model at finite 
temperatures are similar to the ID case. 

For temperatures 0.5J < T < 3J, we derived an ex- 
pression for the mean phonon occupation numbers which 
is of practical value for further investigations of the 
model. 

We investigated the temperature dependence of the 
spin correlation length for two choices of coupling con- 
stants. Our analysis indicates that the model shows long 
range Neel order in the ground state both in the regime 
of weak and strong spin phonon coupling. 
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